## Install & load packages
list.of.packages <- c("maptools","classInt","raster","sp","rgdal","rgeos","PBSmapping","rgl","plot3D","stargazer","mgcv","xtable","mvtnorm"); new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]; if(length(new.packages)){install.packages(new.packages,dependencies=TRUE)}; lapply(list.of.packages, require, character.only = TRUE)


# Custom functions
conSum <- function(x){ifelse(sum(!is.na(x))>0,sum(x,na.rm=T),NA)}
conMean <- function(x){ifelse(sum(!is.na(x))>0,mean(x,na.rm=T),NA)}
conMin <- function(x){ifelse(sum(!is.na(x))>0,min(x,na.rm=T),NA)}
conMax <- function(x){ifelse(sum(!is.na(x))>0,max(x,na.rm=T),NA)}
km2deg <- function(x){x/111.2}
deg2km <- function(x){x*111.2}
deg22km2 <- function(x){x*12365}


######################
## Variable discount rate
######################

# X<-dataX
# d<-1
# k<-1
# r <- .5
# p <- .5
# i <- 2
# polyz <- locations.poly

varyDiscount <- function(X,dvarz,r=NULL,p=NULL,polyz=NULL,t_bat=FALSE){
  
  # Loop over variables
  X.list2 <- lapply(1:length(dvarz),function(d){
    disag <- sort(unique(X$COW.number))
    dvar <- dvarz[d]
    # Loop over wars
    war.list <- lapply(1:length(disag),function(k){#print(paste(dvarz[d],k))
      subz <- X[which(X$COW.number==disag[k]),]
      subz <- subz[order(subz$start_tid),]
      row.names(subz) <- 1:nrow(subz)
      subz.mat <- subz[,c("COW.number","BATTLE_ID","code","initiator")]
      
      if(length(r)>0){
        # current war - discount rate r (over days)
        subz.mat$W_ALL_d <- 0
        if(nrow(subz.mat)>1){
          for(i in 2:nrow(subz.mat)){
            x <- subz[1:(i-1),dvar]
            td <- subz[i,"start_tid"]-subz[1:(i-1),"start_tid"]
            if(t_bat){td <- i-c(1:(i-1))}
            subz.mat$W_ALL_d[i] <- conSum(x/(1+r)^td) 
          }
        }
        
        # current war - same combatant, discount rate r (over days)
        subz.mat$W_SAMEi_d <- 0
        if(nrow(subz.mat)>1){
          for(i in 2:nrow(subz.mat)){
            ix <- intersect(which(subz[,"code"]%in%subz[i,"code"]),1:(i-1))
            x <- subz[ix,dvar]
            td <- subz[i,"start_tid"]-subz[ix,"start_tid"]
            if(t_bat){td <- i-c(1:(i-1))}
            subz.mat$W_SAMEi_d[i] <- sum(x/(1+r)^td,na.rm=T) 
          }
        }
        
        # current war - same opponent, discount rate r (over days)
        subz.mat$W_SAMEj_d <- 0
        if(nrow(subz.mat)>1){
          for(i in 2:nrow(subz.mat)){
            opp.code <- subz[which(subz$BATTLE_ID%in%subz[i,"BATTLE_ID"]&subz$initiator%in%1-subz[i,"initiator"]),"code"]
            opp.bat.id <- subz[subz$code%in%opp.code,"BATTLE_ID"]
            ix <- intersect(which(subz$BATTLE_ID%in%opp.bat.id&(!subz$code%in%opp.code)),1:(i-1))
            x <- subz[ix,dvar]
            td <- subz[i,"start_tid"]-subz[ix,"start_tid"]
            if(t_bat){td <- i-c(1:(i-1))}
            subz.mat$W_SAMEj_d[i] <- sum(x/(1+r)^td,na.rm=T) 
          }
        }
      }
      
      if(length(poly)>0&length(p)>0){
        # current war - same theater
        subz.mat$W_THEATER_d <- 0
        #i <- 3
        for(i in 2:nrow(subz.mat)){#print(i)
          ix <- 1:(i-1)
          bat.id <- subz[ix,"BATTLE_ID"]
          bat.id.i <- subz[i,"BATTLE_ID"]
          bat.id <- bat.id[!bat.id%in%bat.id.i]
          if(length(bat.id)>0){
            bat.id.poly <- locations.poly$BATTLE_ID[which(locations.poly$BATTLE_ID%in%bat.id)]
            subz.poly <- locations.poly[which(locations.poly$BATTLE_ID%in%bat.id),]
            subz.poly.i <- locations.poly[which(locations.poly$BATTLE_ID%in%bat.id.i),]
            wd <- deg2km(gDistance(subz.poly,subz.poly.i,byid=T))+1
            if(length(c(wd))>0){wd.mat <- data.frame(bat.id.poly=bat.id.poly,wd=c(wd))
            wd.mat <- data.frame(bat.id.poly=bat.id.poly,wd=c(wd))
            x <- subz[ix[which(subz[ix,"BATTLE_ID"]%in%bat.id)],dvar]
            x.mat <- data.frame(bat.id=bat.id,x=x)
            x.mat <- merge(x.mat,wd.mat,by.x="bat.id",by.y="bat.id.poly",all=T)
            subz.mat$W_THEATER_d[i] <- sum(x.mat$x*(1/(x.mat$wd^p))/sum(1/(x.mat$wd^p),na.rm=T),na.rm=T)
            }
          }
        }
        
        # current war - same theater, same combatant
        subz.mat$W_THEATERi_d <- 0
        #i <- 3
        for(i in 2:nrow(subz.mat)){#print(i)
          ix <- intersect(which(subz[,"code"]%in%subz[i,"code"]),1:(i-1))
          # x <- subz[ix,dvar]
          bat.id <- subz[ix,"BATTLE_ID"]
          bat.id.i <- subz[i,"BATTLE_ID"]
          bat.id <- bat.id[!bat.id%in%bat.id.i]
          if(length(bat.id)>0){
            bat.id.poly <- locations.poly$BATTLE_ID[which(locations.poly$BATTLE_ID%in%bat.id)]
            subz.poly <- locations.poly[which(locations.poly$BATTLE_ID%in%bat.id),]
            subz.poly.i <- locations.poly[which(locations.poly$BATTLE_ID%in%bat.id.i),]
            wd <- deg2km(gDistance(subz.poly,subz.poly.i,byid=T))+1
            if(length(c(wd))>0){
              wd.mat <- data.frame(bat.id.poly=bat.id.poly,wd=c(wd))
              x <- subz[ix[which(subz[ix,"BATTLE_ID"]%in%bat.id)],dvar]
              x.mat <- data.frame(bat.id=bat.id,x=x)
              x.mat <- merge(x.mat,wd.mat,by.x="bat.id",by.y="bat.id.poly",all=T)
              subz.mat$W_THEATERi_d[i] <- sum(x.mat$x*(1/(x.mat$wd^p))/sum(1/(x.mat$wd^p),na.rm=T),na.rm=T)
            }
          }
        }
        
        # current war - same theater, same opponent
        subz.mat$W_THEATERj_d <- 0
        #i <- 3
        for(i in 2:nrow(subz.mat)){#print(i)
          opp.code <- subz[which(subz$BATTLE_ID%in%subz[i,"BATTLE_ID"]&subz$initiator%in%1-subz[i,"initiator"]),"code"]
          opp.bat.id <- subz[subz$code%in%opp.code,"BATTLE_ID"]
          ix <- intersect(which(subz$BATTLE_ID%in%opp.bat.id&(!subz$code%in%opp.code)),1:(i-1))
          bat.id <- subz[ix,"BATTLE_ID"]
          bat.id.i <- subz[i,"BATTLE_ID"]
          bat.id <- bat.id[!bat.id%in%bat.id.i]
          if(length(bat.id)>0){
            bat.id.poly <- locations.poly$BATTLE_ID[which(locations.poly$BATTLE_ID%in%bat.id)]
            subz.poly <- locations.poly[which(locations.poly$BATTLE_ID%in%bat.id),]
            subz.poly.i <- locations.poly[which(locations.poly$BATTLE_ID%in%bat.id.i),]
            wd <- deg2km(gDistance(subz.poly,subz.poly.i,byid=T))+1
            if(length(c(wd))>0){
              wd.mat <- data.frame(bat.id.poly=bat.id.poly,wd=c(wd))
              x <- subz[ix[which(subz[ix,"BATTLE_ID"]%in%bat.id)],dvar]
              x.mat <- data.frame(bat.id=bat.id,x=x)
              x.mat <- merge(x.mat,wd.mat,by.x="bat.id",by.y="bat.id.poly",all=T)
              subz.mat$W_THEATERj_d[i] <- sum(x.mat$x*(1/(x.mat$wd^p))/sum(1/(x.mat$wd^p),na.rm=T),na.rm=T)
            }
          }
        }
        
      }
      
      subz.mat
    })
    war.mat <- do.call(rbind,war.list)
    war.mat$initiator <- NULL
    names(war.mat)[grep("W_",names(war.mat))] <- paste0(names(war.mat)[grep("W_",names(war.mat))],"_",dvar)
    war.mat
  })
  war.mat2 <- do.call(cbind,X.list2)
  head(war.mat2)
  war.mat2 <- war.mat2[,which(!duplicated(names(war.mat2)))]
  X.w <- merge(X,war.mat2,by=c("COW.number","BATTLE_ID","code"),all.x=T,all.y=F)
  
  return(X.w)
}



######################
## Variable discount rate (inter-war)
######################

# X <- dataX
# dvarz <- "lPOW_max"
# r <- .1
# d <- 1
# i <- 20

varyDiscount_iw <- function(X,dvarz,r=NULL,p=NULL,polyz=NULL,t_bat=FALSE){
  
  # Loop over variables
  X.list2 <- lapply(1:length(dvarz),function(d){
    
    dvar <- dvarz[d]
    subz <- X[order(X$start_tid),]
    row.names(subz) <- 1:nrow(subz)
    subz.mat <- subz[,c("COW.number","BATTLE_ID","code","initiator")]
    
    if(length(r)>0){
      # all wars - discount rate r (over days)
      subz.mat$W_ALL_d <- 0
      for(i in 2:nrow(subz.mat)){
        x <- subz[1:(i-1),dvar]
        td <- subz[i,"start_tid"]-subz[1:(i-1),"start_tid"]
        if(t_bat){td <- i-c(1:(i-1))}
        subz.mat$W_ALL_d[i] <- conSum(x/(1+r)^td) 
      }
      
      # all wars - same combatant, discount rate r (over days)
      subz.mat$W_SAMEi_d <- 0
      for(i in 2:nrow(subz.mat)){
        ix <- intersect(which(subz[,"code"]%in%subz[i,"code"]),1:(i-1))
        x <- subz[ix,dvar]
        td <- subz[i,"start_tid"]-subz[ix,"start_tid"]
        if(t_bat){td <- i-c(1:(i-1))}
        subz.mat$W_SAMEi_d[i] <- sum(x/(1+r)^td,na.rm=T) 
      }
      
      # all wars - same opponent, discount rate r (over days)
      subz.mat$W_SAMEj_d <- 0
      for(i in 2:nrow(subz.mat)){
        opp.code <- subz[which(subz$BATTLE_ID%in%subz[i,"BATTLE_ID"]&subz$initiator%in%1-subz[i,"initiator"]),"code"]
        opp.bat.id <- subz[subz$code%in%opp.code,"BATTLE_ID"]
        ix <- intersect(which(subz$BATTLE_ID%in%opp.bat.id&(!subz$code%in%opp.code)),1:(i-1))
        x <- subz[ix,dvar]
        td <- subz[i,"start_tid"]-subz[ix,"start_tid"]
        if(t_bat){td <- i-c(1:(i-1))}
        subz.mat$W_SAMEj_d[i] <- sum(x/(1+r)^td,na.rm=T) 
      }
    }
    
    if(length(poly)>0&length(p)>0){
      # current war - same theater
      subz.mat$W_THEATER_d <- 0
      subz.mat$W_THEATER_d[1] <- 0
      #i <- 3
      for(i in 2:nrow(subz.mat)){#print(i)
        ix <- 1:(i-1)
        bat.id <- subz[ix,"BATTLE_ID"]
        bat.id.i <- subz[i,"BATTLE_ID"]
        bat.id <- bat.id[!bat.id%in%bat.id.i]
        if(length(bat.id)>0){
          bat.id.poly <- locations.poly$BATTLE_ID[which(locations.poly$BATTLE_ID%in%bat.id)]
          subz.poly <- locations.poly[which(locations.poly$BATTLE_ID%in%bat.id),]
          subz.poly.i <- locations.poly[which(locations.poly$BATTLE_ID%in%bat.id.i),]
          wd <- deg2km(gDistance(subz.poly,subz.poly.i,byid=T))+1
          if(length(c(wd))>0){wd.mat <- data.frame(bat.id.poly=bat.id.poly,wd=c(wd))
          wd.mat <- data.frame(bat.id.poly=bat.id.poly,wd=c(wd))
          x <- subz[ix[which(subz[ix,"BATTLE_ID"]%in%bat.id)],dvar]
          x.mat <- data.frame(bat.id=bat.id,x=x)
          x.mat <- merge(x.mat,wd.mat,by.x="bat.id",by.y="bat.id.poly",all=T)
          subz.mat$W_THEATER_d[i] <- sum(x.mat$x*(1/(x.mat$wd^p))/sum(1/(x.mat$wd^p),na.rm=T),na.rm=T)
          }
        }
      }
      
      # current war - same theater, same combatant
      subz.mat$W_THEATERi_d <- 0
      #i <- 3
      for(i in 2:nrow(subz.mat)){#print(i)
        ix <- intersect(which(subz[,"code"]%in%subz[i,"code"]),1:(i-1))
        # x <- subz[ix,dvar]
        bat.id <- subz[ix,"BATTLE_ID"]
        bat.id.i <- subz[i,"BATTLE_ID"]
        bat.id <- bat.id[!bat.id%in%bat.id.i]
        if(length(bat.id)>0){
          bat.id.poly <- locations.poly$BATTLE_ID[which(locations.poly$BATTLE_ID%in%bat.id)]
          subz.poly <- locations.poly[which(locations.poly$BATTLE_ID%in%bat.id),]
          subz.poly.i <- locations.poly[which(locations.poly$BATTLE_ID%in%bat.id.i),]
          wd <- deg2km(gDistance(subz.poly,subz.poly.i,byid=T))+1
          if(length(c(wd))>0){
            wd.mat <- data.frame(bat.id.poly=bat.id.poly,wd=c(wd))
            x <- subz[ix[which(subz[ix,"BATTLE_ID"]%in%bat.id)],dvar]
            x.mat <- data.frame(bat.id=bat.id,x=x)
            x.mat <- merge(x.mat,wd.mat,by.x="bat.id",by.y="bat.id.poly",all=T)
            subz.mat$W_THEATERi_d[i] <- sum(x.mat$x*(1/(x.mat$wd^p))/sum(1/(x.mat$wd^p),na.rm=T),na.rm=T)
          }
        }
      }
      
      # current war - same theater, same opponent
      subz.mat$W_THEATERj_d <- 0
      #i <- 3
      for(i in 2:nrow(subz.mat)){#print(i)
        opp.code <- subz[which(subz$BATTLE_ID%in%subz[i,"BATTLE_ID"]&subz$initiator%in%1-subz[i,"initiator"]),"code"]
        opp.bat.id <- subz[subz$code%in%opp.code,"BATTLE_ID"]
        ix <- intersect(which(subz$BATTLE_ID%in%opp.bat.id&(!subz$code%in%opp.code)),1:(i-1))
        bat.id <- subz[ix,"BATTLE_ID"]
        bat.id.i <- subz[i,"BATTLE_ID"]
        bat.id <- bat.id[!bat.id%in%bat.id.i]
        if(length(bat.id)>0){
          bat.id.poly <- locations.poly$BATTLE_ID[which(locations.poly$BATTLE_ID%in%bat.id)]
          subz.poly <- locations.poly[which(locations.poly$BATTLE_ID%in%bat.id),]
          subz.poly.i <- locations.poly[which(locations.poly$BATTLE_ID%in%bat.id.i),]
          wd <- deg2km(gDistance(subz.poly,subz.poly.i,byid=T))+1
          if(length(c(wd))>0){
            wd.mat <- data.frame(bat.id.poly=bat.id.poly,wd=c(wd))
            x <- subz[ix[which(subz[ix,"BATTLE_ID"]%in%bat.id)],dvar]
            x.mat <- data.frame(bat.id=bat.id,x=x)
            x.mat <- merge(x.mat,wd.mat,by.x="bat.id",by.y="bat.id.poly",all=T)
            subz.mat$W_THEATERj_d[i] <- sum(x.mat$x*(1/(x.mat$wd^p))/sum(1/(x.mat$wd^p),na.rm=T),na.rm=T)
          }
        }
      }
      
    }
    
    head(subz.mat)
    war.mat <- subz.mat
    war.mat$initiator <- NULL
    names(war.mat)[grep("W_",names(war.mat))] <- paste0(names(war.mat)[grep("W_",names(war.mat))],"_",dvar)
    war.mat
  })
  war.mat2 <- do.call(cbind,X.list2)
  head(war.mat2)
  war.mat2 <- war.mat2[,which(!duplicated(names(war.mat2)))]
  X.w <- merge(X,war.mat2,by=c("COW.number","BATTLE_ID","code"),all.x=T,all.y=F)
  
  return(X.w)
}

#####################################
# set time discount rate (wrapper)
#####################################

# X=dataX
# dvarz=dvarz
# xvarz=xvarz1
# x0=x0
# x1=x1
# runs=runs
# t_bat=FALSE

varyDiscount2 <- function(X,dvarz,xvarz,r=NULL,p=NULL,polyz=NULL,iw=FALSE,t_bat=FALSE){
  
  dvarz1 <- tolower(dvarz)
  xvarz1 <- tolower(xvarz)
  
  # Set discount
  if(!iw){dataX.W <- varyDiscount(X,dvarz,r=r,p=p,polyz=polyz,t_bat=t_bat)}
  if(iw){dataX.W <- varyDiscount_iw(X,dvarz,r=r,p=p,polyz=polyz,t_bat=t_bat)}
  names(dataX.W) <- tolower(names(dataX.W))
  head(dataX.W)
  
  # stardardize coefficients
  st.varz <- c(xvarz1,names(dataX.W)[(grepl("w_",names(dataX.W))&grepl(paste(dvarz1,collapse="|"),names(dataX.W)))|(names(dataX.W)%in%dvarz1[!grepl("ldr_",dvarz1)])])
  temp.mat <- dataX.W[,st.varz]
  for(u in 1:ncol(temp.mat)){temp.mat[,u] <- scale(temp.mat[,u])}
  names(temp.mat) <- paste0("st_",names(temp.mat))
  dataX.W <- cbind(dataX.W,temp.mat); rm(temp.mat)
  #summary(dataX.W)
  return(dataX.W)
}


#####################################
# find optimal discount rate
#####################################
# formz <- as.formula("st_lpow_max~as.factor(cow.number)+as.factor(code)+st_w_samei_d_lpow_max+st_ldeploy_dist+initiator+st_diff_cinc+recruit+st_xpolity+st_opp_xpolity+opp_geneva+start_year+spring+summer+fall")
# formz <- form2b
# 
# X <- dataX
# dvarz <- c("lPOW_max")
# xvarz <- xvarz1
# x0 <- 0
# x1 <- 1

optimDiscount <- function(formz,X,dvarz,xvarz,x0=0,x1=1,runs=3,iw=FALSE,t_bat=FALSE){
  
  # Init
  rz <- seq(x0,x1,length.out = 3)
  rz[1] <- max(rz[1],1e-6)
  
  # loop
  rn <- 0
  while(rn<=runs){
    aicz <- data.frame(r=rz,aic=NA)
    for(r0 in 1:length(rz)){
      dataX.W <- varyDiscount2(X=X,dvarz=dvarz,xvarz=xvarz,r=rz[r0],iw=iw,t_bat=t_bat)
      mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod)
      aicz[r0,"aic"] <- mod$aic
    }  
    aicz
    if(rank(aicz[,"aic"])[1]<rank(aicz[,"aic"])[3]){rz <- seq(rz[1],rz[2],length.out = 3)}
    if(rank(aicz[,"aic"])[1]>rank(aicz[,"aic"])[3]){rz <- seq(rz[2],rz[3],length.out = 3)}
    rn <- rn+1
  }
  return(aicz[which(aicz$aic%in%min(aicz$aic,na.rm=T))[1],"r"])
}


# formz <- form1c
optimDiscount_geo <- function(formz,X,dvarz,xvarz,x0=0,x1=1,runs=3,polyz=NULL,iw=FALSE,t_bat=FALSE){
  
  # Init
  pz <- seq(x0,x1,length.out = 3)
  pz[1] <- max(pz[1],1e-6)
  
  # loop
  pn <- 0
  while(pn<=runs){
    aicz <- data.frame(p=pz,aic=NA)
    for(p0 in 1:length(pz)){
      dataX.W <- varyDiscount2(X=X,dvarz=dvarz,xvarz=xvarz,p=pz[p0],polyz=polyz,iw=iw,t_bat=t_bat)
      mod <- gam(formz,data=dataX.W,family="gaussian"); summary(mod)
      aicz[p0,"aic"] <- mod$aic
    }  
    aicz
    if(rank(aicz[,"aic"])[1]<rank(aicz[,"aic"])[3]){pz <- seq(pz[1],pz[2],length.out = 3)}
    if(rank(aicz[,"aic"])[1]>rank(aicz[,"aic"])[3]){pz <- seq(pz[2],pz[3],length.out = 3)}
    pn <- pn+1
  }
  return(aicz[which(aicz$aic%in%min(aicz$aic,na.rm=T))[1],"p"])
}


######################
## AUC in-sample (GAM)
######################

AUC <- function(y, fitted){
  n1 <- sum(y)
  n <- length(y)
  out <- (mean(rank(fitted)[y == 1]) - (n1 + 1)/2)/(n - n1)
  return(out)}


######################
## AUC out-sample (GAM)
######################

AUC.out <- function(mod,data.out){
  mat.X <- data.out[,names(mod$model)]
  pred.out <- predict.gam(mod,newdata=mat.X[,-1],type="response")
  out <- AUC(y=data.out[,names(mod$model)[1]], fitted=c(pred.out))
  return(out)}

#######################
## RMSE
#######################

rmse <- function(mod){
  sqrt( sum( (mod$y - mod$fitted.values)^2 , na.rm = TRUE ) / nrow(mod$model) )
}

#######################
## GLM predict
#######################

# mod=mod
# var=dfx
# vals=seq(min(dataX[,dfx],na.rm=T),max(dataX[,dfx],na.rm=T),length.out=100)
# fax=list(code=255)

predmanQF <- function(mod,var,vals,quadratic=F,fax=NULL){
  nsim<-10000
  # vars <- unlist(strsplit(as.character(mod$formula)[3],split=" + ",fixed=T))
  vars <- names(coefficients(mod))
  n <- length(vals)
  X <- as.data.frame(matrix(NA,nrow=n,ncol=length(vars)))
  names(X) <- names(coefficients(mod))
  vars. <- vars[!(grepl("as.factor",vars,fixed=T)|grepl("Intercept",vars,fixed=T))]
  X[,vars.] <- as.data.frame(matrix(apply(mod$model[,vars.],2,function(x){median(x,na.rm=T)}),nrow=n,ncol=length(vars.),byrow=T))
  X[,"(Intercept)"] <- 1
  X[,var] <- vals
  if(quadratic==T){
    quadz <- names(mod$model)[grep("I(",names(mod$model),fixed=T)]
    for(k in 1:length(quadz)){
      quadd <- names(mod$model)[apply(as.data.frame(names(mod$model)),1,function(x){grepl(x,quadz[k])})]
      X[,quadz[k]] <- X[,quadd]^2}}
  if(length(fax)>0){
    for(k in 1:length(fax)){
      X[,vars==paste("as.factor(",names(fax[k]),")",fax[[k]],sep="")]<-0
      X[,which(grepl("as.factor",vars)&grepl(names(fax[k]),vars)&(!grepl(fax[[k]],vars)))] <- 0
      X[,which(grepl("as.factor",vars)&grepl(names(fax[k]),vars)&grepl(fax[[k]],vars))] <- 1
    }}
  betas <- rmvnorm(n=nsim,mean=coefficients(mod),sigma=vcov(mod),method="eigen")
  mu <- betas%*%t(X)
  pred.y <- mod$family$linkinv(mu)
  y <- apply(pred.y,2,function(x){median(x,na.rm=T)})
  y.lo <- apply(pred.y,2,function(x){quantile(x,.025,na.rm=T)})
  y.hi <- apply(pred.y,2,function(x){quantile(x,.975,na.rm=T)})
  rrs <- 100*(pred.y[,ncol(pred.y)]/pred.y[,1]-1)
  rr <- mean(rrs)
  rr.lo <- quantile(rrs,.025)
  rr.up <- quantile(rrs,.975)
  return(list(y=y,lo=y.lo,up=y.hi, rr=rr,rr.lo=rr.lo,rr.up=rr.up))}



#######################
## GAM predict
#######################

 # mod=mod
 # var=dfx
 # vals=seq(min(dataX.W[,dfx],na.rm=T),max(dataX.W[,dfx],na.rm=T),length.out=100)
 # fax=list(cow.number=139,code=365)

predmanGAM <- function(mod,var,vals,quadratic=F,fax=NULL){
  nsim<-10000
  vars <- names(mod$model)
  n <- length(vals)
  X <- as.data.frame(matrix(NA,nrow=n,ncol=length(vars)))
  names(X) <- vars
  vars. <- vars[!(grepl("as.factor",vars,fixed=T)|grepl("Intercept",vars,fixed=T))]
  X[,vars.] <- as.data.frame(matrix(apply(mod$model[,vars.],2,function(x){median(x,na.rm=T)}),nrow=n,ncol=length(vars.),byrow=T))
  # X[,"(Intercept)"] <- 1
  X[,var] <- vals
  for(j in 1:length(vars.)){X[,vars.[j]] <- as.numeric(as.character(X[,vars.[j]]))}
  if(quadratic==T){
    quadz <- names(mod$model)[grep("I(",names(mod$model),fixed=T)]
    for(k in 1:length(quadz)){
      quadd <- names(mod$model)[apply(as.data.frame(names(mod$model)),1,function(x){grepl(x,quadz[k])})]
      X[,quadz[k]] <- X[,quadd]^2}}
  if(length(fax)>0){
    for(k in 1:length(fax)){
      X[,which(grepl(names(fax[k]),vars))] <- as.factor(fax[[k]])
    }}
  X <- as.data.frame(X)
  #for(r in 1:ncol(X)){class(X[,names(X[,r])]) <- class(mod$model[,names(X[,r])])}
  names(X) <- gsub(")","",gsub("as.factor(","",names(X),fixed=T))
  preds <- predict.gam(mod,newdata=X,type="link",se.fit=T)
  preds <- cbind(preds$fit,preds$fit-1.96*preds$se.fit,preds$fit+1.96*preds$se.fit)
  preds <- mod$family$linkinv(preds)    
  y <- preds[,1]
  y.lo <- preds[,2]
  y.hi <- preds[,3]
  rr <- 100*(y[length(y)]/y[1]-1)
  rr.lo <- 100*(y.lo[length(y.lo)]/y.lo[1]-1)
  rr.up <- 100*(y.hi[length(y.hi)]/y.hi[1]-1)
  return(list(y=y,lo=y.lo,up=y.hi, rr=rr,rr.lo=rr.lo,rr.up=rr.up))}